
#Set up function to obtain bounds (set concave = 1 if beta_i >= beta_i+1)
cce_lb_ub <- function(concave, weights, acr,lb,ub){
  acr <- abs(acr)
  defaultW <- getOption("warn") 
  options(warn = -1) 
  rm(lprec)
  options(warn = defaultW)
  
  # varCount specifies number of betas
  varCount <- length(weights)
  
  # Initiate LP
  lprec <- make.lp(0,varCount)
  
  # Set objective function as beta_1 + beta_2 ... + beta_n
  objfn <- c(rep(1,varCount))
  set.objfn(lprec,objfn)
  
  # Add constraint to LP that w_1*beta_1 + ... + w_n*beta_n = ACR
  add.constraint(lprec, xt = c(weights), "=", rhs = acr)
  
  # Set lower bound
  lower_bound <- c(rep(lb,varCount))
  upper_bound <- c(rep(ub,varCount))
  set.bounds(lprec, lower = lower_bound, upper = upper_bound, columns = 1:varCount)
  
  # Solve LP if no concavity assumption is made
  if(concave == 0){  
    
    lp.control(lprec,sense='min')
    sol <- solve(lprec)
    # Saves objective function when minimizing objective function
    lb <- get.objective(lprec) 
    # Saves values of betas (note that this isn't unique)
    lb_vars <- get.variables(lprec) 
    
    lp.control(lprec,sense='max')
    sol <- solve(lprec)
    # Saves objective function when maximizing objective function
    ub <- get.objective(lprec)
    # Saves values of betas (note that this isn't unique)
    ub_vars <- get.variables(lprec) 
  
  }
  
  if(concave == 1){ 
    
    concave_Cons <- varCount-1
    
    # Set up concavity constraints
    restrictions <- matrix(0,nrow=concave_Cons,ncol=varCount)
    for (i in 1:concave_Cons){
      restrictions[i,i] <- 1  
      restrictions[i,i+1] <- -1  
    }

    # Add concavity constraints to LP (beta_i - beta_i+1 >= 0) 
    for (i in 1:concave_Cons){
      add.constraint(lprec,restrictions[i,], ">=", rhs = 0)
    }
    
    lp.control(lprec,sense='min')
    sol <- solve(lprec)
    # Saves objective function when minimizing objective function
    lb <- get.objective(lprec)
    # Saves values of betas (note that this isn't unique)
    lb_vars <- get.variables(lprec)
    
    lp.control(lprec,sense='max')
    sol <- solve(lprec)
    # Saves objective function when maximizing objective function
    ub <- get.objective(lprec)
    # Saves values of betas (note that this isn't unique)
    ub_vars <- get.variables(lprec)    
  }
  
  # Output
  result_lp <- list("Lower Bound" = lb, "Upper Bound" = ub, "LB Variables" = lb_vars, "UB Variables" = ub_vars )
  return(result_lp)
}



